#include <math.h>

namespace VGQuantile{

	// P_0 Coefficients (note P_0 is a polynomial of degree 1)
	double c0_0(double lambda, double alpha, double beta)
	{
		return log(1/(alpha - beta))/(alpha - beta);
	}

	double c0_1(double lambda, double alpha, double beta)
	{
		return (-1 + lambda)/(alpha - beta);
	}

	// P_1 Coefficients
	double c1_0(double lambda, double alpha, double beta)
	{
		return ((-1 + lambda)*(-(beta*lambda) + alpha*(2 + lambda) + 2*alpha*log(pow(alpha - beta,-1)))*
		 pow(alpha,-1)*pow(alpha - beta,-2))/2;
	}

	double c1_1(double lambda, double alpha, double beta)
	{
		return pow(alpha - beta,-2)*pow(-1 + lambda,2);
	}

	// P_2 Coefficients
	double c2_0(double lambda, double alpha, double beta)
	{
		return ((-1 + lambda)*pow(alpha,-2)*pow(alpha - beta,-3)*
			(-2*alpha*beta*(-3 + lambda)*lambda + 2*alpha*(alpha*(-4 + lambda) + beta*lambda)*log(pow(alpha - beta,-1)) - 
			lambda*pow(beta,2) + pow(alpha,2)*(-10 + lambda + 2*pow(lambda,2)) - 
			2*pow(alpha,2)*pow(log(pow(alpha - beta,-1)),2)))/4;
	}

	double c2_1(double lambda, double alpha, double beta)
	{
		return ((alpha*(-4 + lambda) + beta*lambda - 2*alpha*log(pow(alpha - beta,-1)))*pow(alpha,-1)*pow(alpha - beta,-3)*
			pow(-1 + lambda,2))/2;
	}

	double c2_2(double lambda, double alpha, double beta)
	{
		return -(pow(alpha - beta,-3)*pow(-1 + lambda,3))/2;
	}

	// P_3 Coefficients
	double c3_0(double lambda, double alpha, double beta)
	{
		return ((-1 + lambda)*pow(alpha,-3)*pow(alpha - beta,-4)*
		 (3*alpha*lambda*pow(beta,2)*(12 + lambda - 3*pow(lambda,2)) + 
		   lambda*pow(beta,3)*(-6 - lambda + pow(lambda,2)) + 
		   3*beta*lambda*pow(alpha,2)*(-46 + 23*lambda + pow(lambda,2)) - 
		   12*alpha*log(pow(alpha - beta,-1))*
			(alpha*beta*(7 - 3*lambda)*lambda - lambda*pow(beta,2) + 
			  pow(alpha,2)*(-14 + 6*lambda + pow(lambda,2))) + 
		   pow(alpha,3)*(188 - 76*lambda - 27*pow(lambda,2) + 5*pow(lambda,3)) - 
		   12*(beta*lambda + alpha*(-5 + 2*lambda))*pow(alpha,2)*
			pow(log(pow(alpha - beta,-1)),2) + 8*pow(alpha,3)*pow(log(pow(alpha - beta,-1)),3)))
		/24;
	}

	double c3_1(double lambda, double alpha, double beta)
	{
		return -(pow(alpha,-2)*pow(alpha - beta,-4)*pow(-1 + lambda,2)*
		  (alpha*beta*(7 - 3*lambda)*lambda + 
			2*alpha*(beta*lambda + alpha*(-5 + 2*lambda))*log(pow(alpha - beta,-1)) - 
			lambda*pow(beta,2) + pow(alpha,2)*(-14 + 6*lambda + pow(lambda,2)) - 
			2*pow(alpha,2)*pow(log(pow(alpha - beta,-1)),2)))/2;
	}

	double c3_2(double lambda, double alpha, double beta)
	{
		return -((beta*lambda + alpha*(-5 + 2*lambda) - 2*alpha*log(pow(alpha - beta,-1)))*pow(alpha,-1)*
		  pow(alpha - beta,-4)*pow(-1 + lambda,3))/2;
	}

	double c3_3(double lambda, double alpha, double beta)
	{
		return (pow(alpha - beta,-4)*pow(-1 + lambda,4))/3;
	}

	// P_4 Coefficients
	double c4_0(double lambda, double alpha, double beta)
	{
		return -((-1 + lambda)*pow(alpha,-4)*pow(alpha - beta,-5)*
		  (3*lambda*pow(beta,4)*(3 + lambda - pow(lambda,2)) - 
			alpha*lambda*pow(beta,3)*(60 + 19*lambda - 26*pow(lambda,2) + pow(lambda,3)) + 
			3*lambda*pow(alpha,2)*pow(beta,2)*(70 + 7*lambda - 35*pow(lambda,2) + 6*pow(lambda,3)) - 
			3*beta*lambda*pow(alpha,3)*(214 - 145*lambda + 2*pow(lambda,2) + 7*pow(lambda,3)) + 
			3*alpha*log(pow(alpha - beta,-1))*(beta*lambda*pow(alpha,2)*(-166 + 109*lambda - 9*pow(lambda,2)) + 
			   lambda*pow(beta,3)*(-6 - lambda + pow(lambda,2)) - 
			   alpha*lambda*pow(beta,2)*(-40 + lambda + 9*pow(lambda,2)) + 
			   pow(alpha,3)*(244 - 156*lambda - 7*pow(lambda,2) + 9*pow(lambda,3))) + 
			pow(alpha,4)*(758 - 471*lambda - 50*pow(lambda,2) + 38*pow(lambda,3) + 4*pow(lambda,4)) + 
			6*pow(alpha,2)*(alpha*beta*lambda*(-23 + 11*lambda) + 3*lambda*pow(beta,2) + 
			   pow(alpha,2)*(52 - 32*lambda + pow(lambda,2)))*pow(log(pow(alpha - beta,-1)),2) - 
			4*(3*beta*lambda + alpha*(-17 + 8*lambda))*pow(alpha,3)*pow(log(pow(alpha - beta,-1)),3) + 
			6*pow(alpha,4)*pow(log(pow(alpha - beta,-1)),4)))/24;
	}

	double c4_1(double lambda, double alpha, double beta)
	{
		return -(pow(alpha,-3)*pow(alpha - beta,-5)*pow(-1 + lambda,2)*
		  (beta*lambda*pow(alpha,2)*(-166 + 109*lambda - 9*pow(lambda,2)) + 
			lambda*pow(beta,3)*(-6 - lambda + pow(lambda,2)) - alpha*lambda*pow(beta,2)*(-40 + lambda + 9*pow(lambda,2)) + 
			4*alpha*log(pow(alpha - beta,-1))*(alpha*beta*lambda*(-23 + 11*lambda) + 3*lambda*pow(beta,2) + 
			   pow(alpha,2)*(52 - 32*lambda + pow(lambda,2))) + 
			pow(alpha,3)*(244 - 156*lambda - 7*pow(lambda,2) + 9*pow(lambda,3)) - 
			4*(3*beta*lambda + alpha*(-17 + 8*lambda))*pow(alpha,2)*pow(log(pow(alpha - beta,-1)),2) + 
			8*pow(alpha,3)*pow(log(pow(alpha - beta,-1)),3)))/8;
	}

	double c4_2(double lambda, double alpha, double beta)
	{
		return -(pow(alpha,-2)*pow(alpha - beta,-5)*pow(-1 + lambda,3)*
		  (alpha*beta*lambda*(-23 + 11*lambda) - 2*alpha*(3*beta*lambda + alpha*(-17 + 8*lambda))*
			 log(pow(alpha - beta,-1)) + 3*lambda*pow(beta,2) + pow(alpha,2)*(52 - 32*lambda + pow(lambda,2)) + 
			6*pow(alpha,2)*pow(log(pow(alpha - beta,-1)),2)))/4;
	}

	double c4_3(double lambda, double alpha, double beta)
	{
		return ((3*beta*lambda + alpha*(-17 + 8*lambda) - 6*alpha*log(pow(alpha - beta,-1)))*pow(alpha,-1)*pow(alpha - beta,-5)*
		 pow(-1 + lambda,4))/6;
	}

	double c4_4(double lambda, double alpha, double beta)
	{
		return -(pow(alpha - beta,-5)*pow(-1 + lambda,5))/4;
	}

	// P_5 Coefficients
	double c5_0(double lambda, double alpha, double beta)
	{
		return ((-1 + lambda)*pow(alpha,-5)*pow(alpha - beta,-6)*
		 (-3*lambda*pow(beta,5)*(60 + 28*lambda - 27*pow(lambda,2) - 2*pow(lambda,3) + pow(lambda,4)) + 
		   30*lambda*pow(alpha,3)*pow(beta,2)*(436 + 24*lambda - 289*pow(lambda,2) + 80*pow(lambda,3) + pow(lambda,4)) + 
		   5*alpha*lambda*pow(beta,4)*(270 + 126*lambda - 143*pow(lambda,2) - 8*pow(lambda,3) + 7*pow(lambda,4)) - 
		   20*lambda*pow(alpha,2)*pow(beta,3)*(249 + 91*lambda - 156*pow(lambda,2) + 6*pow(lambda,3) + 8*pow(lambda,4)) + 
		   5*beta*lambda*pow(alpha,4)*(-6908 + 5796*lambda - 637*pow(lambda,2) - 338*pow(lambda,3) + 35*pow(lambda,4)) - 
		   10*alpha*log(pow(alpha - beta,-1))*(12*lambda*pow(beta,4)*(-3 - lambda + pow(lambda,2)) + 
			  3*lambda*pow(alpha,2)*pow(beta,2)*(-320 + 13*lambda + 148*pow(lambda,2) - 33*pow(lambda,3)) + 
			  alpha*lambda*pow(beta,3)*(258 + 61*lambda - 110*pow(lambda,2) + 7*pow(lambda,3)) + 
			  3*beta*lambda*pow(alpha,3)*(1022 - 855*lambda + 126*pow(lambda,2) + 19*pow(lambda,3)) + 
			  pow(alpha,4)*(-3764 + 3084*lambda - 247*pow(lambda,2) - 200*pow(lambda,3) + 11*pow(lambda,4))) - 
		   pow(alpha,5)*(-37348 + 29842*lambda - 1502*pow(lambda,2) - 2417*pow(lambda,3) + 8*pow(lambda,4) + 
			  77*pow(lambda,5)) + 60*pow(alpha,2)*(beta*lambda*pow(alpha,2)*(-189 + 143*lambda - 20*pow(lambda,2)) + 
			  alpha*lambda*pow(beta,2)*(43 - 4*lambda - 9*pow(lambda,2)) + 
			  lambda*pow(beta,3)*(-6 - lambda + pow(lambda,2)) + 
			  pow(alpha,3)*(296 - 240*lambda + 26*pow(lambda,2) + 8*pow(lambda,3)))*pow(log(pow(alpha - beta,-1)),2) + 
		   40*pow(alpha,3)*(alpha*beta*lambda*(-49 + 25*lambda) + 6*lambda*pow(beta,2) + 
			  pow(alpha,2)*(121 - 89*lambda + 10*pow(lambda,2)))*pow(log(pow(alpha - beta,-1)),3) - 
		   20*(6*beta*lambda + alpha*(-37 + 19*lambda))*pow(alpha,4)*pow(log(pow(alpha - beta,-1)),4) + 
		   48*pow(alpha,5)*pow(log(pow(alpha - beta,-1)),5)))/240;
	}

	double c5_1(double lambda, double alpha, double beta)
	{
		return -(pow(alpha,-4)*pow(alpha - beta,-6)*pow(-1 + lambda,2)*
		  (12*lambda*pow(beta,4)*(-3 - lambda + pow(lambda,2)) + 
			3*lambda*pow(alpha,2)*pow(beta,2)*(-320 + 13*lambda + 148*pow(lambda,2) - 33*pow(lambda,3)) + 
			alpha*lambda*pow(beta,3)*(258 + 61*lambda - 110*pow(lambda,2) + 7*pow(lambda,3)) + 
			3*beta*lambda*pow(alpha,3)*(1022 - 855*lambda + 126*pow(lambda,2) + 19*pow(lambda,3)) - 
			12*alpha*log(pow(alpha - beta,-1))*(beta*lambda*pow(alpha,2)*(-189 + 143*lambda - 20*pow(lambda,2)) + 
			   alpha*lambda*pow(beta,2)*(43 - 4*lambda - 9*pow(lambda,2)) + 
			   lambda*pow(beta,3)*(-6 - lambda + pow(lambda,2)) + 
			   pow(alpha,3)*(296 - 240*lambda + 26*pow(lambda,2) + 8*pow(lambda,3))) + 
			pow(alpha,4)*(-3764 + 3084*lambda - 247*pow(lambda,2) - 200*pow(lambda,3) + 11*pow(lambda,4)) - 
			12*pow(alpha,2)*(alpha*beta*lambda*(-49 + 25*lambda) + 6*lambda*pow(beta,2) + 
			   pow(alpha,2)*(121 - 89*lambda + 10*pow(lambda,2)))*pow(log(pow(alpha - beta,-1)),2) + 
			8*(6*beta*lambda + alpha*(-37 + 19*lambda))*pow(alpha,3)*pow(log(pow(alpha - beta,-1)),3) - 
			24*pow(alpha,4)*pow(log(pow(alpha - beta,-1)),4)))/24;
	}

	double c5_2(double lambda, double alpha, double beta)
	{
		return (pow(alpha,-3)*pow(alpha - beta,-6)*pow(-1 + lambda,3)*
		 (beta*lambda*pow(alpha,2)*(-189 + 143*lambda - 20*pow(lambda,2)) + 
		   alpha*lambda*pow(beta,2)*(43 - 4*lambda - 9*pow(lambda,2)) + lambda*pow(beta,3)*(-6 - lambda + pow(lambda,2)) + 
		   2*alpha*log(pow(alpha - beta,-1))*(alpha*beta*lambda*(-49 + 25*lambda) + 6*lambda*pow(beta,2) + 
			  pow(alpha,2)*(121 - 89*lambda + 10*pow(lambda,2))) + 
		   pow(alpha,3)*(296 - 240*lambda + 26*pow(lambda,2) + 8*pow(lambda,3)) - 
		   2*(6*beta*lambda + alpha*(-37 + 19*lambda))*pow(alpha,2)*pow(log(pow(alpha - beta,-1)),2) + 
		   8*pow(alpha,3)*pow(log(pow(alpha - beta,-1)),3)))/4;
	}

	double c5_3(double lambda, double alpha, double beta)
	{
		return (pow(alpha,-2)*pow(alpha - beta,-6)*pow(-1 + lambda,4)*
		 (alpha*beta*lambda*(-49 + 25*lambda) - 2*alpha*(6*beta*lambda + alpha*(-37 + 19*lambda))*
			log(pow(alpha - beta,-1)) + 6*lambda*pow(beta,2) + pow(alpha,2)*(121 - 89*lambda + 10*pow(lambda,2)) + 
		   12*pow(alpha,2)*pow(log(pow(alpha - beta,-1)),2)))/6;
	}

	double c5_4(double lambda, double alpha, double beta)
	{
		return -((6*beta*lambda + alpha*(-37 + 19*lambda) - 12*alpha*log(pow(alpha - beta,-1)))*pow(alpha,-1)*pow(alpha - beta,-6)*
		  pow(-1 + lambda,5))/12;
	}

	double c5_5(double lambda, double alpha, double beta)
	{
		return (pow(alpha - beta,-6)*pow(-1 + lambda,6))/5;
	}
}
